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The Lagrange-mesh method is a very accurate procedure to compute eigenvalues and eigenfunc- 
tions of a two-body quantum equation. The method requires only the evaluation of the potential 
at some mesh points in the configuration space. It is shown that the eigenfunctions can be easily 
computed in the momentum space by a Fourier transform using the properties of the basis functions. 
Observables in this space can also be easily obtained. 
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I. INTRODUCTION 

The Lagrange-mesh method is a very accurate procedure to compute eigenvalues and eigenfunctions of a two-body 
Schrodinger equation [IHl] as well as a semirelativistic Hamiltonian 0-0]. The trial eigenstates are developed in a 
basis of well chosen functions, the Lagrange functions. Using their special properties, the potential matrix elements 
are simply the values of the potential at mesh points, if they are computed with a Gauss quadrature. At first sight, 
this method could look like a discrete variational method, but this is absolutely not the case since the eigenfunctions 
can be computed at any position. Because of the use of the Gauss quadrature scheme, the method is not variational 
but a great accuracy can nevertheless be reached [Tfj| . The method presented here relies on a mesh of points built with 
the zeros of a Laguerre polynomial, but a general procedure for deriving other Lagrange meshes related to orthogonal 
or non-orthogonal bases has also been developed [Tl|. Even if we only focus on two-body systems in this paper, it is 
worth mentioning that this method can be extended to treat very accurately three-body systems as well in nuclear 
physics as in atomic physics (see for instance Ref. [HI). 

At the beginning, this method was developed in the position space. As we will see below, the potential matrix 
elements are very easy to compute if the interaction is known in terms of the distance r between the interacting 
particles. This is also true for mean values of observables depending on r. For some problems, it can be also useful 
to compute the eigenfunctions in the momentum space by the Fourier transform, as well as observables depending on 
the relative momentum between the particles. We will show that the Lagrange-mesh method can provide these type 
of data very efficiently and very easily, using the fundamental properties of the Lagrange functions. 

The Lagrange-mesh methods in configuration space is described in Sec. [TT1 while Sec. IIIII presents some results in 
momentum space. An ansatz to compute easily the only non-linear parameter of the method is described in Sec. IFVI 
Test calculations are presented in Sec. El and some concluding remarks are given in Sec. I VII 

II. METHOD IN POSITION SPACE 
A. Lagrange functions 

The basic ingredients for the Lagrange-mesh method are a mesh of N points Xi associated with an orthonormal set 
of N indefinitely derivable functions fj(x) [3-0|- The Lagrange function fj(x) satisfies the Lagrange conditions, 

fi{xi) = (1) 

that is to say it vanishes at all mesh points except one. The a;, and Xi are respectively the abscissae and the weights 
of a Gauss quadrature formula 

,•00 N 

/ g(x)dx K,^\ k g(x k ). (2) 
Jo 



* E-mail: gwendolyn.lacroix@umons.ac.be 
t E- mail : claude.semay@umons.ac. be 



Typeset by REVTgX 



2 



As we work with the radial part of wavefunctions, we consider the case of the Gauss-Laguerre quadrature because 
the domain of interest is [0, oo]. The Gauss formula @ is exact when g(x) is a polynomial of degree 2N — 1 at most, 
multiplied by exp(— x). The Lagrange-Laguerre mesh is then based on the zeros of a Laguerre polynomial of degree 
N [l| and the mesh points are given by L jy (xi) = 0. These zeros can be determined with a high precision with usual 
methods to find the roots of a polynomial [13 | (t he Mathematica expression Root does the job efficiently) or as the 
eigenvalues of a particular tridiagonal matrix[l3| . The weights can be computed by the following formula [10] 

N 

InA, = Xi - \nx t + 21nF(7V + 1) - ^ M^; - Xj) 2 . ( 3 ) 

It is worth noting that, for most calculations, it is not necessary to compute the weights A^. The original Lagrange 
functions do not vanish at origin, so it is preferable to use the regularized Lagrange functions whose explicit form is 
given by 

fi{x) = (-l)Xr 1/2 z(z " Xi)- x L N {x) exp(-x/2), (4) 

which is a polynomial of degree N , multiplied by an exponential function. Such a function fi(x) vanishes at the origin 
and at Xj with j =/= i. 

With the Lagrange-mesh method, the solution of a quantum equation reduces (as it is often the case) to the 
determination of eigensolutions of a given matrix. Let us consider the eigenvalue equation 

[T(p 2 ) + V(r)}\i>)=E\iP), (5) 

where T(p 2 ) is the kinetic energy term of the Hamiltonian and V(r) the potential which depends only on the radial 
coordinate r = \ f\. In the following, we will always work in natural units: h = c = 1. A trial state approximation 
of the genuine eigenstate, is expanded on a basis built with these regularized Lagrange functions 

I^>=EQI/,) with (f\fj) = lM^Y lm (r), (6) 
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with f — r/r. The coefficients Cj are linear variational parameters and the scale factor ft, is a non-linear param- 
eter aimed at adjusting the mesh to the domain of physical interest. Contrary to some other mesh methods, the 
wavefunction is also defined between mesh points by (jj) and ©. 

Basis states \fi) built with the regularized Lagrange functions are not exactly orthogonal. But, at the Gauss 
approximation, we have (fj\fi) = So, in the following, all mean values will be performed using the Gauss 

quadrature formula ([2]). In this case, the potential matrix elements are given by 

(fi\V(r)\f j ) = V(hx i )5 ij . (7) 

The potential matrix is both simple to obtain and diagonal. Let us assume that the matrix elements (fi\T\fj) « T^- 
are known. Their computation will be explained in the next section. With ^ and ©, the variational method applied 
to ([5]) provides a system of N mesh equations 

N 

J2 [Tij + V(hxi) % - E fy] C 3 = 0. (8) 
j'=i 

In the Lagrange-mesh method, the Hamiltonian matrix elements are not exactly calculated, but are computed at the 
Gauss approximation. So, the variational character of the method cannot be guaranteed, except if an exact quadrature 
is performed. In practice, for a sufficiently high number of basis states, the method is often variational (eigenvalues 
computed are all upper bounds) or antivariational (eigenvalues computed are all lower bounds) . It has been observed 

1- 3] that the accuracy of the mesh approximation remains close to the accuracy of the original variational calculation 
without the Gauss approximation. So, in most cases, a very high accuracy can be achieved in the framework of the 
Gauss approximation, though the mathematical reasons for the high efficiency of this method are not well known yet 

I"- 

The accuracy of the eigensolutions depends on two parameters: The number of mesh points N and the value of the 
scale parameter h. For a sufficiently high value of N (which can be as low as 20 or 30), the eigenvalues present a large 
plateau as a function of h. This is a great advantage for the Lagrange-mesh method since the non-linear parameter 
must not be determined with a high precision. Nevertheless, if h is too small, a significant part of the wavefunction 
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is not covered by the points of the Lagrange mesh. When h is too large, all points of the mesh are located in the 
asymptotic tail of the wavefunctions and it is then impossible to obtain good eigenvalues. So, it is interesting to have 
a procedure to estimate directly a reasonable value of h in order to avoid a search, which is always time consuming. 
We have remarked that the best results are obtained when the last mesh points are located "not too far" in the 
asymptotic tail. So, if we choose a point r max in the tail of the wavefunction, the value of h can be obtained by 
h = ^max/^AT, where xn is the last mesh point. A procedure to estimate r max will be presented in Sec. IIVI 



B. Kinetic parts 

Let us first look at the matrix P whose elements are P,j = (fi\p 2 \fj)- With $2$, these matrix elements are given by 



h 2 

where I is the orbital angular momentum quantum number, and where 



Uj = I fi{x){~)Mx) dxK-^f'/fa). (10) 



This compact expression is exact for some Lagrange meshes. This is not the case for the regularized Laguerre mesh. 
An exact expression can easily be obtained (see appendix in Ref. 0). However, as shown in Ref. Q, it is preferable 
to use the approximation (|5])-(|10p. The kinetic matrix elements are then even easier to obtain and read [|| 

_ / (-)<-' (.XiXj)- 1 / 2 ^ + Xj )(Xi - x,)- 2 (i £ j), 



Uj - { (l2x 2 )- 1 [4 + {4N + 2)x l -x 2 } {i=j). liiJ 

For a nonrelativistic Hamiltonian, T,j = j^Pij, where \i is the reduced mass of the system. For a more general 

operator T(p 2 ), as the kinetic part of a spinless Salpeter equation 2 y^p 2 + to 2 , the calculation is much more involved. 
The idea is to use a four-step method suggested in Ref. (see also references therein) and applied in Ref. Q: 

1. Computation of the matrix P whose elements are = (fi\p 2 \fj), given by ([5))-(|lip. 

2. Diagonalization of the matrix P. If P is the diagonal matrix formed by the eigenvalues of P, we have 

P = SP D S-\ (12) 
where S is the transformation matrix composed of the normalized eigenvectors. 

3. Computation of T D , a diagonal matrix obtained by taking the function T(x) of all diagonal elements of P D 
(For instance, T(x) — 2\f x + m 2 for the case of a spinless Salpeter equation). 

4. Determination of the kinetic matrix T in the original basis by using the transformation (|12[) 

T = ST D S~ 1 . (13) 

The elements TJy of the matrix computed with (1131) are approximations of the numbers (fi\T(j) )\fj). The calculation 
is not exact for two reasons. First, the elements TSy are computed with an approximate formula ([9|)- (llll) . Second, the 
diagonalization is performed in the limited definition space of the trial function ([5]) . In order to compute exactly the 
matrix elements of the operator T(p 2 ), it is necessary to compute exactly all eigenvalues of the infinite matrix whose 
elements are (T(p 2 )), again exactly computed. This is obviously not possible. It has been shown in Ref. [6|, that this 
four-step procedure can give very good results. 

C. Mean values of radial observables 

The mean value of the operator U(r) for a trial state is given by 

N 

(m(rM) = £ QCjifipir)]^). (14) 
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Using the Lagrange condition (|T|) and the Gauss quadrature ©, this integral reduces to 

N 

(i>\U(rM=J2c]U(h Xj ). (15) 
j=i 

If U is the identity, we recover the normalization condition as expected. A very high accuracy can be obtained with 
this simple procedure p| [l2| . 

III. METHOD IN MOMENTUM SPACE 
A. Fourier transform 

For some particular problems, it can be useful to compute the Fourier transform of a wavefunction in the position 
space in order to obtain the corresponding wavefunction in the momentum space. The Fourier transform 4> FT (p) of 
a wavefunction 4>{r) is defined by 

^ FT (P ) = (2^575 / ^ ) e ^ dft - ( 16 ) 
Using the spherical representation of the wavefunction 

0(f) = Rni(r)Y lm {f), (17) 
and using the spherical expansion of the function e _l ^ r [l6j |. it can be shown that 

<f> FT (p)=R F J( P )Y lm (p), (18) 



where p = \p\ and p — p/p, and where 



R F J(p) = (-l) l \l- I Rni(r)ji(pr)rUr, (!<>) 
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Yi m (p)=i t Y lm (p). (20) 

j; (x) is a spherical Bessel function [n} and Yj m (cc) is called a modified spherical harmonic [l6[ . 
Using expansion the radial part R(r) of the trial function is given by 



The Fourier transform R FT (p) of this radial function is defined by (JTHJ) ■ It is tempting to use the Gauss quadrature 
rule © with the Lagrange condition (TTJ) to perform this calculation. The problem is that spherical Bessel functions 
are rapidly oscillating functions. It is then not obvious that such a procedure could work. Actually, we have checked 
that the Fourier transform of a unique regularized Lagrange function, which is also a rapidly oscillating function, 
cannot be obtained in this way with a good accuracy. Fortunately, the radial part of a wavefunction has a much 
smoother behavior. As we will see on several examples in Sec. El its Fourier transform can be easily obtained in the 
framework of the Lagrange-mesh method by taking benefit of the very special properties of the regularized Lagrange 
function. Using with JT]), the integral (TK)1) simply reduces to 



R* L (p) = (-l)S/-^> Ci y/Xi Xiji(hx iP ), (22) 



where we use the "bar" to indicate that this is not the exact Fourier transform R FT (p). For a sufficiently high value 
of N (which can be as low as 50), R FT (p) Yi m (p) can be a very good approximation of the genuine eigenstate in the 
momentum space for values of p S [0,p max ], where p max can be determined with the procedure used to compute r max 
(see Sec. IIV|) . For values of p > p max , Ji FT (p) can present large unphysical rapid oscillations. These oscillations do 
not develop in R(r), because they are killed by the rapid decreasing of the regularized Lagrange functions. 
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B. Mean values of momentum dependent observables 



The mean value of the operator K(p) for a trial states \tp) is given by 

r°° 9 

(,P\K(PM)= K(P) (R FT (P)) P 2 d P , (23) 
Jo 

where the angular part is already integrated. In this formula, the function R FT (p) can be replaced by R FT (p). 
Good results can sometimes be obtained, but the accuracy cannot be always guaranteed. This is the case when the 
observable grows rapidly with p and needs a very good quality of the asymptotic tail of the wavefunction in the 
momentum space. Actually, it is easier and much more efficient to compute directly 

N 

(mpm = E °i c i (M K (p)\fj)- ( 24 ) 

The matrix elements (fi\K(p)\fj) can be determined by a procedure identical to the one used to compute (fi\T(p 2 )\fj). 
An intermediate step is the calculation of the matrix K , a diagonal matrix obtained by taking the function K(y/x) of 
all diagonal elements of P D (remember that P is linked to the matrix elements of p 2 , notp). The numbers (fi\K(p)\fj) 
are well approximated by the elements of the matrix K obtained by using the transformation <jT2j) : K = S K D S^ 1 . 
As we will see below, a very good accuracy can be reached for the mean values (K(p)). 



IV. SCALE PARAMETER 



An estimation of r max can be computed using the technique developed in Ref. [18j. The first step is to find a 
potential Voo(r) which matches at best the potential V(r) for r — > oo. Three cases are considered in Ref. [ia ]: 

• nr p with k > and p > 0; 

• —n/r p with k > and < p < 1 ; 

• a square well. 

The second step is to choose a trial state |A) which depends on one parameter A, taken as the inverse of a distance. Two 
cases are considered in Ref. [18]: u\(r) oc r l+1 e~ A " r2 / 2 (harmonic oscillator state) and u\(r) oc r l+1 e~ Xr (hydrogen- 
like state), depending on Voair). If the quantum number n is not zero, an effective value of / is used (see Ref. 18]). 
In a third step, the optimal value of A is determined by the usual condition 

^{X\T + V OB {r)\X) =0, (25) 

where T is the kinetic part of the Hamiltonian considered. In the case of complicated T function, the following 
approximation can be used 

(T(p 2 ))^T((p 2 )). (26) 

In particular, we have 

^p 2 + m 2 ^ < y/(p 2 )+m 2 . (27) 

Various expressions for the optimal parameter A are given in Ref. 

Introducing the dimensionless variable s = A r, the regularized radial part u\(s) of the trial state | A) is then analyzed 
to find the value of s e which satisfies the following condition 

UM - e, (28) 



max se[0jOO ] [ux(s)] 



where e (typically in the range 10 4 -10 8 ) is a number small enough to neglect the contribution of u\(s) for values 
of s greater than s e . This is the last step of the procedure, which is very fast and whose details are given in Ref. [l8j]. 
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Note that equation (36) in Ref. [TH has an analytical solution given by (xjy is replaced here by s e in order to match 
the present notations and to avoid a confusion with the last Lagrange-mesh point) 



-(l + l)W- 



(29) 



where W~i is the Lambert function |19j ] and m = 1 or 2 depending on the trial function u\(r). 

At this stage, the ratio s e /X corresponds approximately to a radial distance in the asymptotic tail of an eigenstate 
of the Hamiltonian T + Voo (r) ■ The idea is to identify this distance with the value of r max for the genuine Hamiltonian 
considered. It has been shown in Ref. Q that this procedure works quite well and can give a value of the scale 
parameter h (h = r max /xjv) in the plateau mentioned above. The efficiency of this ansatz is due to the fact that the 
value of h must not be known with a great accuracy in the Lagrange-mesh method. So, a crude determination of r max 
is sufficient and it is not necessary to go beyond the use of the very simple trial functions u\(r) mentioned above and 
the approximation (|26[) for the computation of the kinetic contribution. 

To determine an estimation of p mSLX , let us look at the Fourier transform u^ T (s = p/X) of the trial states considered 
u\(s = Xr): 

u x (s)ocs l+1 e- s2 / 2 => ur(s)oc S l+1 e- s2 / 2 , (30) 
u x ( S )^s l+1 e- s =* ul T (s)cx {s2 S ^ )l+2 - (31) 

If u\(s) is a harmonic oscillator state, m^ t (s) has the same form. So it seems quite natural to set p max = Xs e , since 
both functions present the same ratio (12"51) at the same value of their dimensionless argument. If the trial state is a 
hydrogen-like state, the situation is different since u FT (s) decreases much more faster than u\(s) for large (but not 
too large) values of s. Nevertheless, the simple choice p max = Xs e works quite well also, as it will be shown below. 
So, finally, we have 

r ma x = Se/A and p ma x = As e , (32) 
with s e and A determined by the procedure described above. 

V. NUMERICAL TESTS 

In this section, several tests will be performed for the Lagrange-mesh method with both nonrelativistic and semirel- 
ativistic kinematics. We will focus on the quality of wavefunctions and observables in the momentum space since the 
efficiency of the method in the position space has already been demonstrated elsewhere In order to estimate 

more precisely the quality of the Fourier transform (|22|). we define a "quality factor" Q(p*) 



Q(p*) — max 
pe[0,p»] 



(p)-u^(p) 



m aXpG[o,p»] I mFT (p)I 

where u FT (p)/p = R^J(p) given by (|22|) and u FT (p)/p = R FT (p) is the exact solution in momentum space. 

A. Confining semirelativistic Hamiltonian 

Let us consider the ultrarelativistic two-body system with a quadratic potential 

H = 2Vp2 +ar 2 . (34) 

This Hamiltonian is particularly interesting because it is probably the only one with a semirelativistic kinematics 
which is (partly) analytically solvable. With an appropriate change of variable, this Hamiltonian can be recast into 
the form of a nonrelativistic Hamiltonian with a linear interaction [20j, for which solutions are known for S-states. 
The eigenvalues for / = are given by 

E n0 = (4a) 1 / 3 |a„|, (35) 
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where a n is the (n + l)th zero of the Airy function Ai [l7|. The corresponding regularized eigenfunctions are obtained 
directly in the momentum space |21] 



u no(p) =pRno(p) 




Ai (a„) 



Let us note that Ai 2 (s) ds — Ai' (a n ). Using the generalized virial theorem [22], it can be shown 
i' 2 |nO) = (nO\ar 2 \nO) where |n0) is a S-eigenstate. Moreover, all powers of p can be computed exactly [23| 



(nO\-\/p' 
we have: 



(36) 

that 
. So, 



(37) 

(nO\p 4 \nO) = (|) 4/3 J|(8|a n | 4 + 25|«„|). (38) 

To perform the following calculations, we have set a = 0.25. The units of the results are given in powers of the unit 
chosen for the only energy scale of the system a 1 / 3 . Using the Lagrange-mesh method with N — 10 and e = 10~ 4 , 
the eigenvalues (|3"5j) can already be obtained with a relative error smaller than 1%. But, to obtain a good Fourier 
transform of the wavefunction, it is necessary to use more points. As we can see on Fig. [TJ the agreement can be 
very good for the main part of u FT (p). With N = 20, unphysical oscillations appear just before p max . With N = 40, 
they develop halfway between p max and 2p max . With N — 80 (not presented here), the asymptotic behavior is correct 
till 2p max . In these 3 cases, for which e = 10~ 8 , we have respectively Q(p max ) — 0.034, 0.0042, 0.0052. The quality 
factor first decreases rapidly due to the improvement of the wavefunction for large values of p, and then stabilizes 
because the quality of the wavefunction stays constant in the low-p part. It is possible to improve the quality factor 
by decreasing the value of e (increasing the value of p max ). For N — 40, the value of Q(p ma *) decreases from 0.015 to 
0.0020 when e varies from 10~ 4 to 10~ 12 . 



«"(/>) 
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FIG. 1: The exact solution (|36|l with a = 0.25 for n = 2 is compared with the corresponding approximation given by formula 
(|22l) for < p < 2p max . The value of p max is determined with the procedure presented in Sec. IIVI with e — 10 -8 . 



Some observables for a particular eigenstate, I = and n = 2, computed with formula ([M]) are presented in Table U 
and compared with the exact values. Similar results are obtained for other eigenstates. A very good accuracy can 
be obtained with a quite small number of points. Actually, it appears that the precision does not automatically 
increases with N. On the contrary, for a given value of e, the accuracy is optimal for a given number of points. This 
behavior is typical of semirelativistic Hamiltonians. This is due to the computation of the kinetic part which requires 

Our experience is that an 



a supplementary approximation than the use of the Gauss quadrature rule (see Sec. Ill B 
optimal value for an observable can be found by looking at extrema or plateau in the behavior of this observable as 
a function of N for a given value of e. In the next section, we will see on an example that accuracy increases with N 
for a nonrclativistic system. 



B. Hydrogen atom Hamiltonian 



We consider now a completely different case, the hydrogen atom: the kinematics is nonrelativistic and the Coulomb 
potential, —ct/r, is non-confining. The eigensolutions in the position space are well known and their Fourier transform 
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TABLE I: Some observables with a — 0.25 for the eigenstate I = and n = 2, computed with formula (1241) and compared with 
the exact values. Results are given in powers of the unit for a 1 / 3 . 



Exact 



{VP) jgj (exp(-pVa 2 / 3 )} 
1.84019 (a) 24.0273 (6) 0.109740 (c) 



e = 10" 6 AT = 



10 1.84198 
20 1.84265 
40 1.84399 



23.6260 
24.0735 
24.0982 



0.108562 
0.109299 
0.108892 



10" 8 N 



10 1.81901 
20 1.84163 
40 1.84236 



23.6006 
24.0545 
24.0680 



0.112181 
0.109512 
0.109359 



(a) 



Computed with l|37|l: { ' Computed with 



(c) 



Computed with quadrature using (|36 



can be expressed in term of the Appell Hypergeometric function F% [24[ . As these special functions are difficult and 
lengthy to obtain accurately, it is more convenient to work with numerically c omp uted eigensolutions in momentum 
space. Particular momentum dependent observables can be exactly computed |23| : 



(P 2 ) = 

(p 4 > = V 



(n + l + 1) 2 ' 
4 8n + 21 + 5 
' (2; + l)(n + Z + l) 4 ' 



(39) 
(40) 



where r\ = [ia, with fi the reduced mass. 

To perform the following calculations, we have set mi — 940 MeV, — 511 KeV, a = 1/137. The units of the 
results are given in powers of keV. Some observables for a particular eigenstate, / = 1 and n = 1, computed with 
formula (|24[) are presented in Table [TT] and compared with the exact values. Similar results are obtained for other 
eigenstates. Again, a very good accuracy can be obtained with a quite small number of points. This time, accuracy 
always increases with N for a given value of e, as already found in previous studies 



sr otpc 
[2, 10], 



TABLE II: Some observables for the hydrogen atom eigenstate / — 1 and n — 1, computed with formula (1241) and compared 
with the exact values. Results are given in powers of keV. 



Exact i.54414 (a) 11.9218 (f,) 0.786997 (c 

e = 10~ 6 N = 10 1.54417 11.9225 0.787043 
20 1.54414 11.9218 0.786995 
40 1.54414 11.9218 0.786994 

e = 10~ 8 N = 10 1.54711 11.9471 0.787255 
20 1.54414 11.9218 0.786997 
40 1.54414 11.9218 0.786997 

^ Computed with ((39j) ; < - b - 1 Computed with ()40[): ( - c - ) Computed with quadrature of the numerical Fourier transform of the 

wavefunction in position space. 



A good Fourier transform of the main part of the wavefunction u FT (p) can be obtained with a small number of 
points, around N = 20-40. But, to obtain a good asymptotic tail, it is necessary to use more points, as we can see on 
Fig. [2] With N — 100, unphysical oscillations appear before p max . With N = 200, they develop halfway between p max 
and 2p max . For e = 10~ 6 , we have respectively Q(p max ) = 0.504, 0.097, 0.00028, for N = 50, 100, 200. Nevertheless, 
the quality factor Q(p») can be as small as 10 -6 if is in the main part of the wavefunction. It is also possible to 
improve the quality factor by decreasing the value of e (increasing the value of p max ) ■ 
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FIG. 2: The accurate numerically computed (exact) Fourier transform of the hydrogen atom wavefunction for I = 1 and n = 1 
is compared with the corresponding approximation given by formula (|22[1 for < p < 2p max - The value of p max is determined 
with the procedure presented in Sec. IIVI with e = 10~ 6 . 



VI. CONCLUDING REMARKS 

The Lagrange-mesh method is a procedure to compute eigenvalues and eigenfunctions of quantum equations. It 
is very simple to implement and can yield very accurate results for a lot of observables, specially for nonrelativistic 
kinematics. At the origin, the method has been developed in the position space since the evaluation of potential 
matrix elements requires only the computation of the interaction at some mesh points. This is due to the use of 
a Gauss quadrature rule with the fact that the basis functions satisfy the Lagrange conditions, that is to say they 
vanish at all mesh points except one. Using this very special property, we have shown that the computation of the 
wavefunction in the momentum space by the Fourier transform of the wavefunction in the position space can be easily 
performed with a very good accuracy. Moreover, mean values of momentum dependent operators can also be easily 
and accurately calculated using a technique similar to the one used to compute the semirelativistic kinetic matrix 
elements. This shows again the great efficiency of the Lagrange-mesh method which can yield very accurate results 
for a minimal computational effort. We can wonder if this technique could also be used directly in the momentum 
space, for instance in the case where the interaction is only known as a function of the relative momentum. This 
question will be addressed in a subsequent paper. 
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